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ABSTRACT 


The purpose of this research is to design, simulate and 
implement a robust autopilot system for the vertical mode of 
operation of the Archytas prototype. Archytas is an Unmanned 
Air Vehicle that is designed to take off and land vertically, 
and to transition to horizontal forward flight. A feedback 
control scheme is designed for both the single-input, single- 
output and the multi-input, multi-output subsystems using 
optimal control techniques. In this research, the linear 
quadratic regulator performance measure is modified to allow 
for its application to the tracking problem solution. 
Additionally, the control systems are designed using reduced 
order models. Computer simulations show that the reduced 
order controller designs provide results comparable to the 
full order controller designs. Successful hardware tests with 
roll rate control system validated the reduced order model 


design philosophy used in this research. 
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I. INTRODUCTION 


A. THE ARCHYTAS CONCEPT 

The current inventory of Unmanned Air Vehicles (UAVs) is 
not able to meet the expanding need for real-time intelligence 
at the Marine expeditionary unit or Army battalion level. The 
Naval Postgraduate School UAV Flight Research Lab is directing 
efforts at developing a ducted-fan vertical takeoff and 
landing (VTOL) vehicle to meet these increasing needs. The 
NPS air vehicle, named Archytas, is serving as a technology 
demonstrator to evaluate the concept of a winged ducted-fan 
VTOL aircraft. The research is being directed at applying the 
technology and equipment developed in the U.S. Marine Corps’ 
Airborne Remotely Operated Device (AROD) program and the U.S. 
Army’s AQUILA progran. 

Archytas, pictured in Figure 1.1, is designed to take off 
and land vertically. After climbing to altitude, Archytas 
(irene) tien tomhorizomtal flight by pitching about its 
center of gravity to a wings level attitude. The positioning 
of the duct and wings (including the canard) allow for the 
vertical takeoff and landing capability. The ability to 
transition to horizontal flight will extend the vehicle’s 


range and provide the capability for a high speed dash. 
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Figure 1.1 Sketch of Archytas 


B. THE CONTROL PROBLEM 

Archytas is powered by a. vertically mounted 28-horsepower 
engine turning a three-bladed propeller. The use of a single 
propeller in a duct (ducted fan) simplifies the design, but 
intensifies the stability and control problems. The dynamic 
behavior about a given axis is coupled with other vehicle 
dynamics. In particular, there are three types of coupling: 


1. The single propeller design introduces a gyroscopic 
coupling between the pitch and yaw axes. 


2. Reactive torques are applied to the roll axis as the 
engine speed is varied. 


3. A loss of lift due to thrust occurs when the vehicle is 
pitched during the translation to horizontal flight. 


It was the goal of this thesis to design a control system 
which would allow stable flight of Archytas during takeoff and 
landing, and during the vertical mode of operation. Because 
of the coupled nature of the Archytas control problem, linear 


quadratic regulator control theory was used. 


C. THESIS ORGANIZATION 

In Chapter II, the general nonlinear equations of motion 
are developed. These equations of motion are then applied to 
Archytas with special attention given to the effects of the 
spinning propeller. Additionally, the equations of motion are 
linearized about the hover operating condition using the 


small-disturbance theory. 


Chapter III discusses the state space representation of a 
general system and develops a procedure for selecting state 
variables for the optimal control tracking problem. Key to 
the design of the Archytas control system is the proper 
selection of the state variables. 

In Chapter IV, three control law designs are formulated 
based on the linearized hover equations. One design is for 
the single input, single output (SISO) roll rate controller. 
The second design is for the SISO altitude rate controller. 
These two SISO systems are similar in their development. 
Next, the multiple input, multiple output (MIMO) pitch and yaw 
angle controller is designed. Central to each control law 
design is the use of a reduced order model to simplify the 
design process and physical implementation. Finally, these 
control laws are applied to the nonlinear model for 
validation. 

Chapter V discusses the field test results of thewimeus 
rate controller. The roll rate controller was evaluated with 
Archytas mounted on a test stand to allow a roll about the 
longitudinal axis. In addition, conclusions based on the 
computer simulations, the field tests, and recommendations for 


future research are presented. 


II. MODELING THE ARCHYTAS PROTOTYPE 
The purpose of this chapter is to develop a suitable 
dynamic model of the Archytas prototype. Because Archytas 
is a ducted fan device, special attention must be given to 


the significant gyroscopic contribution of its propeller. 


A. DERIVATION OF RIGID BODY EQUATIONS OF MOTION 

The rigid body equations of motion in this section are 
developed for the Archytas prototype in the following way. 
First, Archytas is regarded as a single rigid body, and the 
equations of motion are derived with respect to a set of 
body fixed axes. These equations are the general equations 
governing aerodynamic flight for all aircraft. Next, the 
changes introduced by the spinning rotor are evaluated and 
included in the equations. These are the gyroscopic effects 
due to the propeller. Finally, the development of a 
complete model is undertaken for the specific case of 
Archytas using the actual measurements and experimental data 
from the AROD prototype as first approximations. 

1. Force and Moment Equations 

The general equations of motion are developed for a 

typical aircraft in References 1 and 2. A combination of 
the two approaches is taken here to arrive at the set of 


equations describing Archytas. The equations of motion are 


obtained from Newton’s second law, which states: The 
Summation of all external forces acting on a body is equal 
to the time rate of change of the momentum of the body; and 
the summation of the external moments acting on the body is 
equal to the time rate of change of the moment of the 
momentum (angular momentum). The time rates of change of 
linear and angular momentum are referred to an absolute or 
inertial reference frame. This absolute or inertial 
reference frame is an axis system fixed to the Earth. 
Figure 2.1 depicts both the body fixed axes and the inertial 
reference frame.[(Ref. 1: p.84] 

Newton’s second law can be expressed in the 


following vector equations: 


(2.1) 


= es (2.2) 
~ d 


where F is the externally applied force, M the externally 
applied moment about the center of mass, v the velocity 
vector, and H the angular momentum vector about the center 
of mass. 

The vector equations, in scalar form, consist of 
three force equations and three moment equations. The force 


equations can be expressed as 


"ee 





wee eZ IFZF#@ng, . 
aes 


y a 


Xf 


Fixed Frame 


Y¢ 


Figure | 


Figure 2.1 Body and Fixed Axes System 


aaa (mw) ; (2.3) 


«at ae 20 Gils 


where F,, F,, F,,and u, v, w are the components of the fomee 


y? 
and velocity along the x, y and z axes, respectively. The 
force components are composed of contributions due to the 
aerodynamic, propulsive, and gravitational forces acting on 
the aircraft, the Archytas prototype for the purpose of this 


thesis. The moment equations can be expressed in a similar 


manner: 


d a d (2.4) 


where L, M,N ancdeae, Jol wemeol, 20 c= the components of the moment 
and moment of momentum along the x, y and Z axes 
respectively. 

Now, considering Figure 2.2, let 5m be an element of 
mass of the Archytas prototype, v be the velocity of the 
elemental mass relative to the inertial axes, and let 6F be 
the resultant force that acts upon it. Newton’s second law 


then gives the equation of motion of dm: 
dy 
6 F=5m— . (2.5) 
= dt 


The total force acting on the vehicle is a summation of all 


the forces that act upon all the elements. The internal 





Figure 2.2 An Element of Mass on the Archytas 


forces, those exerted by one element upon another, all oceimie 
in equal and opposite pairs, by Newton’s third law, and 
hence contribute nothing to the summation. Thus LOF=F is 
the resultant external force acting upon the vehicle. The 


velocity of the differential mass 6m is given by 


(2.6) 


ale 


vey 
C 


were v, is the velocity of the center of mass of the 
aircraft and dr/dt is the velocity of the element relative 
to the center of mass. Substituting this expression for the 


velocity into Newton’s second law, equation (2.1), yields 


ar 
DoreF=y mS (y +—) . (2.7) 
BED tn ee gE) 
Assuming that the mass of the aircraft is constant, equation 


(2.7) can be rewritten as 





_ @ ar 2.8 
Eads (Yet qe) om (Te 
OE: 
a 2 
agp = ré6m. (2.9) 
ia aE eee 


Because r is measured from the center of mass, the summation 


Xrém is equal to zero and the force equation (2.9) becomes 


10 


dv 
ay 


(2.10) 
lowes 


=m 





which relates the external force on the aircraft to the 
motion of the vehicle’s center of mass. 

The relationship between the external moment and the 
rotation of the aircraft is obtained from a consideration of 
the moment of momentum. For the differential element of 


mass, 6m, the moment of momentum is by definition 5H=rxvém. 


The moment equation can be written as 
6M=— 6H=—= (rxv) dm. (2.11) 


The velocity of the mass element can be expressed in terms 
of the velocity of the center of mass and the velocity of 


the mass element relative to the center of mass: 


VV +——=V +@xl ; (2.12) 


where w is the angular velocity of the vehicle and r is the 
position of the mass element measured from the center of 


mass. The total moment of momentum can be written as 
EP Vey (Zxv_) dm+)" [rx (wxzr)]6m . (2yai'3 ) 


The velocity v, is a constant with respect to the summation 


and can be taken outside of the summation Sign 


aioe 


H=()° rdm) x vty (rxwxz) 6m . (2.14) 


The first term in equation (2.4) YS zero because Che Cemm 


Xrbm=0 as explained previously. The position vectors and 
angular velocity can be expressed as 


Taxi + yaaa; (2.15) 


@=pitqitrK ; (2.16) 


where p, q and r are the scalar components of w, and i, j, k 
are unit vectors in the directions of x, y and Zz. 
Substituting w and r into equation (2.14) and expanding, H 


can be written as 
H=(pitqitrk) }> (x?+y?+z?) &m-)” (xit+tyit+zk) (px+qytrz) 6m 


(25270 
The scalar components of H are 
Hi= DY (¥7+27) pO ec ain ae) eee 


H,=-pyxy dm+qy(x?+z?) 6m -ryyz 6m ; (2.18) 


H=-p) x2 6m-qQ) yo Omtz) (ety orn 


The summations in the above equations are the mass moment 


and products of inertia of the aircraft and are defined as 


2 


reffftyrretiom, tof ffsy om 
Figzalali\nec?+.22).dm Iaz={ { [xz Om ; (20829 ) 
mene sty") 5m, 1y,-{{fyz Om. 


The terms I,, I, and I, are the mass moments of inertia of 
the body about the x, y and z axes, respectively. The terms 
with the mixed indices are called the products of inertia. 
Both the moments and products of inertia depend on the shape 
of the body and the manner in which its mass is distributed. 
The larger the moments of inertia the greater the resistance 
the body will have to rotation. Applying the notation of 
equation (2.19) to equation (2.18), yields the scalar 


equations for the moment of momentum: 


pee ee 7 


Hy lees) le. Ley (2.20) 
fe eeee an 


If the reference frame is not fixed to the aircraft, then, 
as the aircraft rotates, the moments and products of inertia 
will vary with time. To avoid this difficulty an axis 
system will be fixed to the aircraft (body axis system). 

Now the derivatives of the vectors v and H referred to the 


rotating body frame of reference must be determined. 


a8 


It can be shown that the derivatives of an arbitrary 
vector A referred to a rotating body frame having an angular 
velocity w can be represented by the following vector 


identity 








where the subscripts I and B refer to the inertial and body 
fixed frames of reference, respectively. Applying this 


identity to equations (2.1) and (2.2) yields 








dy 

F=m—<| +m(wxyv_) ; (2.22) 
clita c 
we +W xH . (2.23) 
=" dials 


These are the general equations governing aerodynamic flight 


and have the scalar components: 


Fi=m(u+qw-rv), F,=m(v+ru-pw), F,=m(w+pv-qu); (2.24) 


L=H,+QH,-rH,,  M=H,+rH,-pH,, N=H,+pH,-qH, . (2am 
The components of the force and moment acting on the 
aircraft are composed of aerodynamic, gravitational and 
propulsive contributions. 

At this point, it is recognized that most aircraft 
have a plane of symmetry. If the xz plane is selected to 
coincide with this plane of symmetry, then from equation 
(2.19), I,=I,=0 must be satisfied. However, for the case of 


xy yz 


14 


the Archytas, the three blades of the propeller provide two 
planes of symmetry, xz and xy. Thus, the products of 


0 


inertia, I . and I,,, equal zero. The moment equations 


xy / 
for Archytas can now be written as 

Baemaaia() — i) % 

Mone mt I) ; (2.26) 


Bree oc 1 .— 1.) 


2. Effect of the Spinning Rotor 

Archytas, like AROD, is a gyroscope. The single 
propeller rotates about the longitudinal vehicle axis to 
produce a downwash or jet of air through the duct which 
makes up the Archytas body. This spinning rotor exerts a 
gyroscopic moment on the body of the vehicle. Reference 3 
states that in developing the equations of motion for 
aircraft with propellers which exert gyroscopic moments on 
the body "more often than not, such gyroscopic moments turn 
out to be negligible." However, as demonstrated by Bassett 
[Ref. 4: p. 19], in the AROD case the angular momentum of 
the propeller, Ijapzop, equals 11.3 ft’-lb,/sec. Compared 
with AROD’s nominal total mass of 2.64 lb, (85 lb,), it 1s 
clear that the angular momentum imparted by the propeller is 
significant and that gyroscopic effects will play a large 
part in modeling the dynamic behavior of AROD. For 
Archytas, I, is identically equal to I,agop, 11.3 ft*-1b,/sec. 


Compared with Archytas’ nominal total mass of 3.11 1b, (100 


BES) 


ib,), it is clear that similar to AROD thewgyroscopic 
effects will be significant and must be included in modeling 
the dynamic behavior of Archytas. This gyroscopic moment 
can be accounted for as follows. 


Angular momentum, H, due to the propeller is 


—pl 


defined as 
(2.27) 


where I, is the propeller moment of inertia, and W, 153m 
propeller angular velocity. Since the propeller lies in the 
yZ plane and spins symmetrically about the x axis, H, is 


directed only along x and ithe = Ons Equation (2.27) becomes 


ies Dato ati (2.28) 
Etkin (Ref. 2: p.93] states that the resultant angular 
momentum of an aircraft with spinning propellers is obtained 
Simply by adding H, to the H previously defined iy equation 
(2.20). Adding equations (2.28) and (2.20) and keeping in 


mind that the products of inertia equal zero, yields 


Has bi Loe 
Hy= cds (2.29) 


7 a a 


Applying equation (2.29) to equation (2.25), the moment 
equations can be written for the specific case of Archytas 


as 
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JS Ossi ys as Maal ek 
Meee i 15) trl 0, j (2.30) 


Nee OG l= 2,.).- GI 2Os : 


and the force equations are those of equation (2.24). 
3. Orientation and Position 

Because the frame of reference developed for the 
equations of motion is fixed to the aircraft, and moves with 
it, the position of the aircraft cannot be described 
relative to it. The orientation and position of the 
aircraft can be defined in terms of a fixed frame of 
reference as shown in Figure 2.3. [Ref. 1: p. 89] 

The orientation of the aircraft can be described by 
a series of three consecutive rotations, whose order is 
important. The angular rotations are called the EF£uler 
angles. The orientation of the body frame with respect to 
the fixed frame can be determined in the following manner. 
The aircraft is imagined first to be oriented so that its 
axes are parallel to the fixed frame and the following 
rotations are then applied. [Ref. 2: pp. 89-91] 


1. A rotation WV about oz,, carrying the axes to Cx,yZ, 
(Bringiamg Gx to its fimal azimuth) . 


2. A rotation @ about oy,, carrying the axes to Cx;y32; 
(bringing Cx to its final elevation). 


3. A rotation $ about ox,, carrying the axes to their 


final position Cxyz (giving the final angle of bank to 
the wings). 


po) 


X4,X 
Fr in -: Ww 
=) e- 4 Pe) _ Y2:¥3 
Oe ~. wae Cy & ®D q 
—— eS 
= 2 
i e 
/ 4 2 
ae 
_ on 
® «. | wet 
Ze | Flight path 
8 Z1,Z2 y; 
r x 
O 
Earth-fixed axes 
Vertical 
Ze 


Figure 2.3 Relationship between Body and 
Fixed Axes System 
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Now that the Euler angles are defined, the flight velocity 
components relative to the earth-fixed reference frame can 
be determined. To do this, let the velocity components 
along the x,, y,;, z, frame be dx/dt, dy/dt, dz/dt and, 
Similarly, let the subscripts 2 and 3 denote the components 
along x,, Y), Z, and x3, y;, Z,, respectively. From Figure 


pee, itwean be shown that 


dy dz 
———— =] Z — =V, ; —— —W. 4 (2.31) 
a cl Cima” 

where 
u=u,cosY—v_siny ; 
V,=u,Ssin+v,cosyp ; (2.32) 
W,=W, ; 

and 
u,=u,cosO+w,sind ; 
oa emr (2.33) 
W,=-u,SinO@+w,cos6 ; 

and 


U,=U ; 
v,=vcos®-sin® ; (2.34) 


w,=vsin®+wcos® ; 


from this, the absolute velocity in terms of the Euler 


angles and velocity components in the body frame can be 
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determined. Note the shorthand notation S,=sinV, C,=cosY, 


Sg=sin8, etc, used in the following equations: 


EX) (OsCy SeSiCi Cog VOASeS arse =e ie 


dt 

d 

ae =| C5Sy “SpcpSyt Coee ene yoann (2.35) 
az 

| de| | ~%e 59Ce Colo = | LY 


Integration of these equations yields the aircraft’s 
position relative to the fixed frame of reference. The 
relationship between the angular velocities in the body 
frame (p, gq and r) and the Euler rates ( 9, w and ®) can 


also be determined from Figure 2.3. 


Gao’ Ce CsS5 Q (2.36) 


Equations (2.36) can be solved for the Euler rates in terms 
of the body angular velocities and is given by equation 


(22367) 


| [1 stand Cytand] [> 
O)=|9 Ce “So | |q (29am 
| |OnS—gs€eO CysecO) |r 


By integrating the above equations, the Euler angles (0, ¥ 


and ®) can be determined. 
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4. Gravitational and Thrust Forces 

The gravitational force acting on the aircraft acts 
through the center of gravity of the aircraft. Because the 
body axis system is fixed to the center of gravity, the 
gravitational force will not produce any moments. However, 
the gravitational force will contribute to the external 
force acting on the aircraft and will have components along 
the respective body axes. From Figure 2.4 the gravitational 


force component in the direction of each axis is found to be 


X,=-mg cosOcos® ; 
Y,=- mg cos@Osin? ; (2.38) 


Z,=-mg sind . 


With the aerodynamic forces (including the propulsive 
forces) denoted by (X, Y, Z), the resultant external forces 


are 


F,=X-mg cos@Ocos¥; 
A= Yemgieosusin® ; (2.39) 


F=Z-mg sind . 


S.- Summary of Equations of Motion 
In the previous sections, the equations that 
completely describe the dynamic behavior of Archytas have 
been developed. Equations (2.24) and (2.30) define the 
externally applied forces and moments which are represented 


by F,, F,, F, and L, M, N. Through the Euler angles, the 


yl 
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9 





Figure 2.4 Components of Gravity acting 
along Body Axes 


behavior of Archytas can be observed relative to the Earth. 
Specifically, the translational velocities for the fixed 
frame of reference, dx/dt, dy/dt and daz/dt, can be 
determined from the body-fixed velocities, u,v,w and the 
Euler angles, 9, VW and $, using the transformation of 
equation (2.35). Additionally, equation (2.37) describes 
the relationship between the Euler angles and the body 
angular velocities, p, q and r. Table 2.1 gives a summary 


of the rigid body equations of motion. 


B. ARCHYTAS NONLINEAR SYSTEM EQUATIONS 

1. Applied Forces and Moments 

The Archytas model is developed using the 

measurements and data from the AROD prototype as first 
approximations. These measurements and data were obtained 
by the AROD project engineers based on wind tunnel tests. 
This data consists of tabular results that describe the 
aerodynamic lift and drag coefficients and physical 
measurements of constants such as weight, moments of inertia 
and servo gains. This tabulated data forms the basis from 
which the applied forces and moments may be determined. The 
forces and moments, which are computed from the data are of 
two types: aerodynamic and thrust. This data is listed in 


Appendix A. 
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TABLE 2.1 SUMMARY OF GENERAL EQUATIONS OF MOTION 


X-mg CyCy = m(U+qw-rv) 

Y+mg CgSe = m(v+ru-pw) Force equations 
Z-mg Se = m(w+pv-qu) 

ESi DYqr (ie 1, aor 

M= Tq+op (1, a gear Moment equations 


N=1 7 +*pgt lt, -2) ae eae 
p= -¥S, 
G=OCy+ POS Angular velocities 


T=? Coen -OS, 


0=QCy-ISs 
P=p+qSgTo+ICeTe Euler rates 


P= (gSy+rCy) secd 


Velocity of aircraft in the fixed frame in terms of Euler 


angles and body velocity components 


EX! (CyCy SgSpCy-CoSy CySpCy* SoS] cy 


=|CeSy SeSgSytCaCy CoSeSy-SeCy Vv 


ae 
ay 
dt 
dz 
dt 


— Se SeCe Cole 





a. Total Angle of Attack and Body Roll Angle 
The aerodynamic data describes the forces and 
moments relative to the vehicle’s total velocity vector, 
Vrotr) These forces and moments (in the body-fixed coordinate 
system) depend on the total angle of attack (a’) and the 
body roll angle (A). The total angle of attack and body 


roll angle can be defined in terms of the velocity 
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components as shown in Figure 2.5. The equations for a’ and 


A are given as: 


« “=sint Weer © (2.40) 
u 
and 
A=tan7t” ' (2.41) 
W 


b. Aerodynamic Forces and Moments 
The forces computed from the tabular data are 
lift (F,) and drag (F,). These forces are depicted in Figure 
2.6(a). The transformation from lift and drag to F,,, F, and 


F, 1S given as 














F.x| | sina “ Soa Py (2.42) 
F.y.| @iCOsa “ -sime Jr ,| 
and 
Pay a sinA (2°43) 
mee icosA| ° 











where F,,, F,, and F,, are the forces in the x, y and z 
directions due to the aerodynamic data. 

Similarly, the aerodynamic moments applied to the 
body axes as a result of the vehicle’s movement through the 
air can be derived. These moments (shown in Figure 2.6(b)) 
are referred to as the aerodynamic angular moments of roll 


(L,), pitch (M,) and yaw (N,) and are given as 
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y, -Z, plane 





, 





wo) 
N 


Angle of Attack (a’)and 
Body Roll Angle (A) 


Figure 2.5 
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iene scoso, | 5, 
F,,| {[-cosa’ sina’ || F, 





(a) Forces 





(b) Moments 


Figure 2.6 Aerodynamic Forces and Moments 


Pa 

















Le | | sing -cosa “| |M, (2.44) 
My2| |-COSsa —Simler jae 
and 
Ms _ {sind (2.45) 
Nel |\COS Aye 











where M,, and M, are the yaw and roll moments relative to 
Vrotr> My, 18 the moment in the y-z plane. The relationships 
that define Ff), Fy wand M,, M, were developed by the AROD 
project engineers and are listed in Appendix A. [Ref. 15] 
c. Control Forces ar Moments Due to Command Inputs 

The forces and moments previously discussed are a 
result of the vehicles motion through the air. A second 
category includes those forces and moments which are a 
result of the commanded inputs. These inputs control the 
rotor speed and the displacement of the control surfaces. 
[Refs 4:5 pen 

(1) Forces Due to Induced Thrust of the Ducted 
Fan. The commanded inputs include the ability to change the 
rotor speed; thus, changing the thrust provided by the 
ducted fan. Due to the orientation of the body-fixed axes, 
the force due to the thrust (F,,) is directed completely 


along the x-axis. The relationship that defines F,, was 


developed by the AROD project engineers and is given as 
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F,-=THOVER + XRPM * 6... (2.46) 


where THOVER is a constant equal to the nominal thrust at 
hover. THOVER is set equal to Archytas prototype weight of 
76.5 ft/lbs. XRPM is the siope of the thrust versus engine 
rpm curve. For the hover rpm of 712.0943 rad/sec (6800 


rpm), XRPM is equal to 0.2387 lb,/rad/sec. 46-.,, is the 


change in engine rpm from the nominal hover rpm in rad/sec. 
(2) Moment Due to Ducted Fan Effects. A 
gyroscope imparts no torque on its axis if it spins with a 
constant angular rotation. In hover (constant rotor speed), 
Archytas behaves similarly to a gyroscope. If the rotor 
accelerates (positively or negatively) a torque is applied 
to the axis. This torque is accounted for in Equation 
(2.30). However, Archytas is a ducted fan and the drag 
between the rotor tip and the inside body wall creates a 
moment about the x-axis (roll). The project engineers for 
AROD determined an approximation for this moment based on 
experimentation. Because the Archytas duct is identical to 
the AROD duct, the moment determined for AROD applies to 


Archytas and 1S given as 


Lon = ian rpm (2.47) 
where 6,,, is defined above. Ky, is a constant which 1s 


dependent on the duct geometry and is equal to 0.0729. (Ref. 


4: pp. 33-34} 


ao 


(3) Moments Due to Control Surface Displacements. 

The commanded inputs also include the displacement of the 
four control surfaces (vanes) within the downwash from the 
duct. The displacement of these vanes within the downwash 
of the duct imparts moments about the body axes. Figure 2.7 
Shows how the vanes are symmetrically arranged below the 
duct. The vanes are displaced by a servo mechanism 
connected directly to the top of each vane. Vanes (1) and 
(3) are operated together as "elevators" and impart a moment 
about the y-axis (pitch). Vanes (2) and (4) together are 
the "rudder" and contribute a moment about the z-axis (yaw). 
Vanes (1) and (3) displaced in opposite directions and (2) 
and (4) displaced oppositely work as "ailerons" to impart a 
moment about the x-axis (roll). The actual torque applied 
by each combination of vanes was determined experimentally 
and described by "constants of effectiveness" which were 
calculated by the AROD project engineers. These constants 
of effectiveness can be applied directly to Archytas. [Ref. 
4: p.34] 

The constants of effectiveness are given the symbols 
Lice Meeep ANG N,, for their contribution of moments about the 
roll, pitch and yaw axes due to a displacement by the 
ailerons, elevator and rudder. The relationships resulting 


in moments about the three body axes are 
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P, 6) (2.48) 


DaGlOr 


as Se e / 


Neale VpseearO rf 


~ 


where 6,, 6, and 6, are the displacements of the aileron, 


elevator and rudder respectively. Pro, and Yiu, are scaling 
factors. The relationships that define P,,.,, and Yr, were 
developed by the AROD project engineers and are listed in 
Appendix A.  Lyae, Meee and N.- equal -150,379.57, -112,716.87 
and -128,774.80, respectively, with units of l1b,-in’- 
rad/sec*. Because of the symmetry of Archytas, cross 
coupling of the control surfaces is negligible and is 
lgnored. {Ref. 4: pp. 35-36] 
2. Servo Equations for Control Surfaces and Throttle 
The model airplane servos used to drive the control 
vanes and throttle linkage are identical and can be modeled 
aS second-order dynamical systems. The response of these 
Servos to a step response was measured. These measurements 
were used to compute the natural frequency and damping 
ratio. The results are summarized below in the form of 
constants H, and H,. A detailed explanation as to how these 
results were obtained is contained in Appendix B. 
a. Control Surface Servos 
Three equations will describe the operation of 
the elevators, rudders and ailerons. For each of these 


equations, at least two servos are operating at the same 
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time on different control vanes. Each servo receives a 
command input (pulse-width modulated signal) which results 
in an angular displacement of the servo. The corresponding 


differential equations for the servos can be written as 


Ome 0 29,6 +H, u; / 
C= 240 -H,O0+H,u, ; (2.49) 
6 .=-H,6 pd, 0 US / 


where H, and H, are the servo gain constants equal to 71.1 


oeme745.8. u,, U. and u. are the servo inputs. _6,, 46, and 
6. are the servo position angles. 


b. Throttle Servos and Engine 
The servo motor used to open and close the 
throttle is identical to the servo motors used for the 
control surfaces. The 2-cycle, 2-cylinder gasoline engine 
with dual carburetors can be modeled as a first order lag 


system. The complete third order system is 


Dees tO) eta.) anaes 
b rpm £9 pom? O KS 7 ; 
where w, 1S the lag time constant equal to 2.0 rad/sec and 
Kz 18 a scaling factor equal to 837.758 rad/sec/rad. 
3. Summary 
The result of the previous section was nine 


equations completely describing the dynamic behavior of 


Archytas. These equations, combined with the equations that 
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define the servos and engine, form the complete nonlinear 
system. This nonlinear system will be the base model 
throughout this thesis. Table 2.2 lists the equations 
rearranged so that the dynamic variables of interest may be 
determined. Additionally, the applied forces and moments 
defined in this section have been substituted into the 
equations. These applied forces and moments complete the 


development of the model for the specific case of Archytas. 


C. LINEARIZING THE ARCHYTAS MODEL 

There exist many analytical and graphical techniques for 
controller design and analysis of linear systems. 
Conversely, there are no good methods available for solving 
a wide class of nonlinear systems. Thus, in the design of 
control systems it is practical to first design the 
controller based on the linear system model generated by 
neglecting the nonlinearities of the system. These 
nonlinearities are neglected by linearizing the model about 
a steady-state reference condition. The designed controller 
is then applied to the nonlinear system model for validation 
and subsequent redesign if necessary. In this section, a 
linear model is generated from the equations of Table 2.2 
based on steady-state assumptions and physical 
approximations. (Reis ss pices) 

The nonlinearities of the Archytas system equations fall 


into two categories: (1) nonlinear combination of states 
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Table 2.2 ARCHYTAS NONLINEAR SYSTEM EQUATIONS 


u=(rv-qw) -gcosécos¥ += Pie (Pomel) tel re (9 vam ld 


V= (pw-ru) +gcosOsin¥+—F,,(F;, Fy) 


w= (qu-pv) -gsinO+=—F,, (aa ee) 


8=gcos®-rsin® 
6 =p+gsin®tand+rcos®tané 
YP =(gsin®+rcos®) secé 


(e.g. Equation (2.30)) and (2) discontinuous functions (e.g. 
table lookup of aerodynamic force and moment coefficients). 
These nonlinearities will be neglected and the nonlinear 
equations will be replaced with linear approximations. 
1. Steady-State Assumptions 
The motion of Archytas in the hover mode consists of 


small perturbations from a steady-state condition. The 


3 


steady-state hover condition is defined as a’ equal to 90°. 
All translational and angular movement is very small. The 
steady-state hover condition results in the following 
simplifications of the nonlinear system equations. 


1. The product of two small numbers is an extremely small 
number, thus terms involving the products of 
translational or angular velocities are equal to zero 
(e.g. rq, pq, (pw-ru) are set equal to zero). 


2. The aerodynamic forces (F,,, Fy, Fy) and moments (L,, M,, 
N,) are very small for a’ equal to 90°, and can be 
neglected in the hover flight condition. (Note: 
Because the vanes lie within the downwash of the 
propeller, the vane effectiveness coefficients can not 
be neglected.) 


3. The sine of a state is equal to the state and the 
cosine of a state is equal to one. This is the small 
angle approximation for angles less than 15 degrees. 


4. K,,, is a small number equal to 0.0729. For hover or 
very small translational velocities, Oe is a very 
Small number. Therefore, the product of Ky, and O Fae 
is neglected. 


5. The propeller angular velocity, ,, is considered 
constant. Thus, the propeller angular acceleration, 


®,, 1s equal to zero. 
Table 2.3 lists the Archytas system equations when the above 
simplifications are applied to the nonlinear system 
equations. Note that much of the coupling between states 
and all of the nonlinear products of states have been 


eliminated. 
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Table 2.3 ARCHYTAS LINEARIZED HOVER EQUATIONS 


} 7-H, S .-H,5 ,.+H,uU, 


6 rpm —@ £9 pom? @ pK pb e 





2. Physical Approximations of Force and Moments 
Memwnonec ie ine taole 2.3 1s a function of the 
engine rpm and is computed by equation (2.46). The moment 
terms, L, M, and N,, are functions of the displacement of the 
alleron, elevator and rudder. They are computed by equation 
I 


(2.48). The moments of inertia, I iad L_, were 


xf yl 


determined by the AROD project engineers and are listed in 


Appendix A. The angular velocity of the propeller, W,, in 


ET, 


the hover mode equals 712.0943 rad/sec (6800 rev/min). The 
weight of the AROD prototype is equal to 76.5 lbs. Using 
these constant parameters, the forces and moments for the 
hover steady state reference condition were computed and are 
listed in Table 2.4. Substituting the results of Table 2.4 
into the equations of Table 2.3 yields the linearized hover 


equations, which are summarized in Table 2.5. 


D. SUMMARY 

The nonlinear differential equations of motion of a 
rigid body were developed from Newton’s second law. The 
equations were then modified for the specific case of 
Archytas. This modification included the effect of the 
spinning three-bladed rotor. The equations were linearized 
about a steady-state hover condition. Next, the applied 
forces and moments, and the servo and engine equations were 
defined. The forces and moments specific for the hover 
reference condition were computed and applied to the 


Archytas linear model. 
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Table 2.4 FORCES, MOMENTS AND CONSTANTS 


Horces 
=71 8 Fmnpy? THOVER + XRPM 8-55 


aeff“a Sat aed 


Lp=-150,379.578, Fm, -=76.5+0.23878.25 


M. efter er! fantor® 2 


M,=(-112,716.87) (1)8 


e 


ae er! factor Or 


N,=(-128,774.80) (1)6 


Eq 


Mements of Inertia 
i -7062. 3908, 1-7 %6e.22, 71-7729 .58,e1,,=60.552 (ib, in*) 


Angular velocity 


W)=712. 0943 222. 


ec 


We] Vi m 
welght=76.5 lbs. ; gravity=32.174 = ; 
sec 


Ena; Time ¢ eee Neumaai 


=a 0 =f, K,=837.758 rad/sec/rad 
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Table 2.5 LINEARIZED HOVER STATE EQUATIONS 






= 150,379.57 5 - 53 298, 
7063.39 


q=~--—__ [-112,716.876,-52, 440, 822) =-14 vse, seen 


TENG eee 
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sii tA OPTIMAL CONTROL THEORY 
The purpose of this chapter is to develop a procedure 
for applying optimal control theory to the solution of 
tracking problems. First, several reasons for desiring an 
optimal solution are presented. Next, the state space 
representation (both continuous and discrete) of a system is 
given. Finally, the application of optimal control to the 


solution of tracking problems is illustrated. 


A. WHY OPTIMUM CONTROL FOR ARCHYTAS? 

The first reason for seeking an optimum controller is 
that feedback gains can be computed for a much broader range 
of control problems. Specifically, optimal control provides 
solutions for high order, multiple-input, multiple-output 
(MIMO) systems. Such systems are often intractable with 
classical methods. The pitch and yaw angle controller for 
Archytas is an eight order MIMO system. Thus, optimal 
control is the preferred method. 

Additionally, optimal control lends itself nicely to a 
discrete time solution of the control problem. Archytas 
will employ an on board digital computer to perform inflight 
Stability and control. While a continuous time controller 
can be easily discretized in many cases, design of a sampled 


data controller will simplify the procedure.[Ref. 3: p. 58] 
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Finally, optimal control is known to provide robust and 
insensitive solutions to the feedback control problem. 
Assuming that an appropriate performance measure is chosen 
to determine the optimal feedback gain matrix, K, the 
solution can be expected to have a fair degree of tolerance 
to plant model inaccuracies. Clearly, robustness 1s not 
only a desired property of the controller, it is an absolute 
necessity if the controller is to be applicable to both the 


linear and nonlinear models of Archytas.{Ref. 3: Dp. 58-59} 


B. STATE SPACE REPRESENTATION 

The state of a system may be defined as the minimum 
amount of information required such that (given the input to 
the system) the response of the system is completely 
determined for all future time. For dynamic systems, the 
response of the system is defined by the differential 
equations that model the system, the initial conditions and 
the forcing function. The number of state variables or 
states is equal to the total order of the systems 
differential equations. In order to provide a systematic 
mathematical approach to analysis of the system, it is 
convenient to describe the system by a set of first-order 
differential equations. This set of equations is called the 
state equations and constitute the basis for the state space 


representation. [Ref. 8: pp. 206-207] 


42 


1. Continuous Time Systems 
The state space representation of a general n” order 
continuous-time, time-invariant system is described by the 


following matrix state equations: 


ieeeen xc )ems UG) ; (3.1) 


WARE =O (6). (3.2) 


where the definitions in Table 3.1 apply to a system with 2 
control inputs and m measurable outputs. The equations 


TABLE 3.1 
STATE SPACE DEFINITIONS FOR CONTINUOUS-TIME SYSTEMS 


Dimension Definition 
State vector 
Output vector 
Plant matrix 


Gentzred distribwutien matrix 


GuEpUut Disteloutmen mati 





listed in Table 2.5 are linear and time-invariant. 
2. Discrete Time Systems 
As was noted earlier, there are many benefits for 
seeking a discrete time solution for the Archytas control 
problem. Therefore, the automatic control systems designed 
will focus on the application of optimal control theory to 


discrete time systems. 
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Similar to the continuous-time system, the state 
Space representation of a general n” order discrete-time 


system is described by the following matrix state equations: 


X(K+1) =® x(k) +E u(k) ; (3.3) 
VK) =C sail (3.4) 
® and [ are defined as: 
@=e* ; (3.5) 
[=f ‘e4'de B ; (3.6) 


where T is the sampling period and k is an integer time 
index. Reference 10 provides a detailed development of the 
relationship between continuous-time and discrete-time 
systems. The definition in Table 3.2 apply to a system with 
£ control inputs and m measurable outputs. 


TABLE 3.2 
STATE SPACE DEFINITIONS FOR DISCRETE-TIME SYSTEMS 


Dimension Definition 
oo) State vector 
ol") Output vector (0 <~m<n) 
x mn) Plant matrix 


x 0) Control distr2 buia onumaiteaae: 


Ine 
* 


(k) 
y (k) 
® 
iP 
g 


awe) Output Distribution matrix 
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cs LINEAR QUADRATIC REGULATOR PROCEDURE 

The theory behind the linear quadratic regulator is well 
developed. Procedures exist which make the design of 
controllers using linear quadratic regulator theory easily 
obtainable. The purpose of this section is to illustrate 
the application of the linear quadratic regulator to the 
tracking problems associated with Archytas. Reference 11 
provides a detailed development of the linear quadratic 
regulator theory. 

des Quadratic Cost Function 

The discrete LOR synthesis problem is that of 

determining the control that minimizes the performance 


MeaSure ; 


J=S° x7(k) QO x(k) + uT(k) Ru(k) ; (3.7) 
k=0 


where 
a = Scalar cost of operating the system; 


x(k) = State vector at discrete times; 

u(k) = Control vector at discrete times; 

Q = Symmetric state trajectory weighting 
Maree 1x ° 

R = Symmetric control weighting matrix; 

T = Matrix transpose operator. 
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The solution to this optimization probleme tie ened 


controller: 
u(K) =-(RB* Myth) = Ae (3578) 
where M satisfies the algebraic Riccati equation (ARE): 


0 = Ae - eRe (3.9) 


Ite 


(Ref. 11: pp. 345-346] Many software packages, including 
the program MATLAB used in this thesis, contain subroutines 
that calculate the value of K for a given dynamic system and 
performance measure. 
2. Performance Weighting Matrices 

The optimal control is fully specified by the system 
model and the weighting matrices Q and R. Q penalizes 
deviation of the state vector x from the origin and R 
penalizes the use of too much control effort. Trial and 
error was used in selecting values for these performance 
weighting matrices. The relationship between the weighting 
matrices Q and R and the dynamic behavior of the closed-loop 
system depends of course on the matrices A and B and is 
quite complex. The approach taken in the design of the 
controllers for Archytas was to solve for the gain matrices 
K that result from a range of weighting matrices Q and R, 
and then simulate the corresponding closed-loop response. 
The gain matrix K that produced the response closest to 


the desired design objectives was chosen. [Ref. 11: p. 341] 
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3. Optimal Tracking Systems 

The regulator and tracking problems appear very 
Similar. The tracking problem attempts to drive the states 
of the system to a desired level; whereas, the regulator 
attempts to drive all of the states to zero. Their 
differences present conceptual difficulties and sensitivity 
problems when viewed in a practical context [Ref. 12: pp. 
643-647]. Burl (Ref. 13] presents a comprehensive 
development of the subtleties encountered when applying LOR 
synthesis to the tracking problem. This development is 
generalized below to demonstrate these subtleties and to 
indicate a procedure for selecting the proper form of the 
performance measure. The following development is applied 
to a first order system, but the results are applicable to 
systems of arbitrary order. 


Given the scalar plant: 


ya ayn Gjerou we). (3.10) 
The purpose of the control system is to drive the output 
y(t) to the constant reference input r. This results in the 
error equation: 
Cie eV (if) (3512) 


where the desire is to drive e(t) to zero. The application 
of linear regulator theory requires that this error be a 


linear combination of the states of the system. This can be 
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readily accomplished by generating a new state equation for 


e(t): 
eit) —==V(E) =ay (2) - Saeeae (3.12) 
adding and subtracting Ar yields: 


é(t).=Ay(.t)+4An=ome Buea 
é(t) =A(-rty(t)]-Bu(t)+Ar ; 
é(t)=-Ae(t)-Bu(t)+Ar . (3) 


Since the error equation, Equation (3.19), is linear and 


time invariant, the performance measure: 


J=[(0e?(t) tu? (t)}de ; (3.14) 
0 


should provide an optimal solution for the system of 
Equation (3.21). However, as shown by Burl [Ref. 13] this 
optimal control problem will not yield a solution. To 
demonstrate this fact, note that the existence of this 
integral requires that both u(t) and e(t) approach zero as t 
tends to infinity. If e(t) approaches zero, then e(t) must 


also approach zero which from Equation (3.21) yields: 


0=-A(0) -Bu(t) tar = u(t)=Sr 4 0 (3 Ss) 


This nonzero value for u(t) implies that the performance 
measure is infinite. A solution to this problem would be to 


formulate the performance meaSure as 
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oo 


Zz 
J=f [oe?(t} +{u(t) -4r) ldt . (3.16) 
9 


This will result in the existence of a theoretically optimal 
control provided that the model of the system (the 
parameters A and B) is exactly known. When errors exist in 
the model, a nonzero steady state error will exist. The 
resulting control system will be mathematically correct, but 
it will be unacceptable for many applications due to its 
Sensitivity to changes in the plant parameters. The control 
system will be super-tuned (it will not be robust).[{Ref. 12: 
p. 645] 

This undesirable result can be overcome by letting the 
controller find the appropriate steady state value of the 
control. This is achieved by application of the performance 


measure: 


J=[ (Qe? (t) +u?(t)] de ) (3.17) 
0 


The control that minimizes Equation (3.25) given the system 
of Equation (3.21) is found by first differentiating 


Equ@derome(3.21), yielding: 
Sie Acie ea Buiit). . (3.18) 


Noting that the error is the integral of é(t) results in 


the state space system: 
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2 ares (3.19) 
=) 


+ 

















This system and the performance measure of Equation (3.25) 
form a linear regulator problem whose solution is state 


feedback: 


(Cc) 
ysis 


(3.20) 





u(t) =-[k, Kile 


The gains k, and k, are computed by application of linear 
regulator theory. The actual control that is applied to the 


system is found by integrating Equation (3.28): 


5 
u(t) =-k,fe(t) dt-k,e(e) (3.21) 


0 


The resulting controller incorporates proportional plus 
integral feedback of the error. This fact is reasonable 
since the system to be controlled, Equation (3.21), 1s of 
type 0 with a steady state disturbance input. 

To summarize, when LOR synthesis is applied to the 
tracking problem, the proper choice of state variables will 
help to eliminate sensitivity problems and ensure system 
robustness. The state variables should not contain any 
input dependent terms, and they should asymptotically 
approach zero. The system tracking error and its time 
derivative are natural state variable for the LQR synthesis 


procedure. In essence, the key to a successful design 
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depends upon the proper formulation of the performance 
measure. The performance measure er Equation’ (3.25) 
represents one possible formulation that yields reliable 
results. This performance measure is used in the design of 
the Archytas control systems. These ideas are extended to a 
general multiple input, multiple output (MIMO) system by 


Burl (Ref. 13}. 


D. SUMMARY 

In this chapter, the application of the linear quadratic 
regulator to the tracking problem has been addressed. A 
procedure for formulating a performance measure applicable 
to tracking problems was developed. Additionally, the state 
Space representations for both continuous-time and 
discrete-time systems were presented. In the next chapter, 
the results of this chapter are applied to generate control 


systems for Archytas. 


oa 


IV. CONTROL SYSTEM DESIGN FOR ARCHYTAS 
The purpose of this chapter is to use optimal control 
theory to design an automatic flight control system for 
Archytas during the hover mode of flight. Because of the 
coupled nature of the Archytas control problem, a linear 


quadratic regulator (LOR) synthesis technique is pursued. 


A. ARCHYTAS CONTROL SUBSYSTEMS 
The automatic control system is logically separated into 
three subsystems according to the linearized equations of 
Table 2.5. The three control subsystems are: 
1. Roll rate controeiien:: 
2. Altitude rate controller; 


3. Pitch angle and yaw angle controller. 


Because each of these control subsystems is designed 
independently, any cross-coupling which may occur between 


the subsystems is not taken into account. 


B. ARCHYTAS ROLL RATE CONTROLLER 
1. The Roll System 
The roll rate controller for Archytas will serve 


three primary functions: 


Sa 


1. Allow the operator to command a desired rotation 
velocity about the vehicle’s longitudinal, or xX, axis. 
This rotation velocity will permit the operator to 
position Archytas by terminating the roiling motion 
when a desired angle is achieved. This capability is 
criticai during the landing phase in order to position 
Archytas correctly with respect to the wind. 


2. Eliminate the rotation velocity imparted to Archytas 
from the effect of cross-winds. 


3. Eliminate the rotation velocity imparted to Archytas 
from the reactive torques applied to the roll axis as 
the engine speed is varied. 


The simplified equation of motion which describes 


Archytas’ roll rate subsystem is given as: 
p=-21.298, . (4.1) 


The aileron servo dynamics are modeled in Appendix B as a 


second order dynamical system with a natural frequency, ,, 
of 52.4 rad/sec (8.34 Hz) and a damping ratio, ¢, of 0.707. 
The corresponding differential equation is: 

Of ero, - 2745, 8042745. eu... (4.2) 


From Equations (4.1) and (4.2), the state equation 


can be written in the matrix form X=A x+B u: 


P| {9 -21.29 0 |P QO.) 
o 3) =o 0 1 |.) + 0 u, (ome )) 
6. 0 -2745.8 -71.1 8. Ah va sy. 3S) 


The roll system tracking error is defined as: 
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B6=D. ee, (4.4) 


where p. 1s the input command and p is the measured roll 
rate. From Equation (4.4), the differential equation for 


the tracking error is: 


E=-p. (4.5) 


If Equation (4.5) is combined with the time derivatives of 
Equations (4.1) and (4.2) a new state equation can be formed 


that is appropriate for tracking system design (as discussed 


in Chapter III). This state equation is: 
Ep} fo -1 0 0 y*p 0 
Py j0 0-21.29 0 ee (4.6) 
$.} |0 0 0 1 44, 0 2 
§| [0 0 -2745.8 -74.1}§ | [2745.8 
L a! 


2. Roll Rate Controller Design 
a. Sampling Frequency Selection 
The application of optimal control theory to 
discrete-time systems requires that an appropriate sampling 
frequency be determined. The sampling frequency, f,, is 
Simply the inverse of the sampling period, T. A general 
rule of thumb in control systems is to sample a system such 


that 
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16a ,0 aan (4.7) 


where f,,, 1S the maximum bandwidth of the system. From 
Appendix B, the natural frequency of the servos is 52.4 
rad/sec (8.34 Hz). Bassett [Ref. 3: p. 76] demonstrates how 
the Archytas MIMO system can be considered as a set of 
several SISO systems. Each system has a different bandwidth 
for the purpose of determining a sampling frequency. The 
subsequent highest natural frequency is equal to 4.64 
rad/sec (0.74 Hz). Because the natural frequency of the 
servos is a factor of ten greater than the highest natural 
frequency of 4.64 rad/sec, the servo natural frequency will 
be used to determine the required sampling period, T. From 
Equation (4.7), £,=10(52.4 rad/sec)=524.0 rad/sec (83.4 Hz), 


using 83.4 Hz as the sampling rate yields 


J=— =——— =0.012 sec . 4.8 
ie "Corre ( 


For the controller designs of this thesis, a sampling period 
of 0.01 seconds will be used. 
b. Discretizing the Roll Rate System 
MATLAB provides tools for computing the discrete- 
time matrix state equation. With the sampling period T 
equal to 0.01 seconds, the discrete-time state space can he 


written in the matrix form x(k+1)=® x(k) +E u(k): 
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ED (<1) 11-0 caunenon ‘ieee to | 0 
DVR a) Gee 9.2048 -9.0008/ P(x) | |-0.00e1]. 
= + u(k). (4.9) 
Fae coe ae (oe 0.8935 0.0067 16,(k)| | 0.1065] ° 
8.(k+1)}/ 10 0 -18.5258 013936 |§o(ay{ [18-5258 


c. Gain Determination 
The optimal control is determined from the state 


equation and the performance measure: 


2 pis) 
J=S))[Bp(K) B(K) 8,(K) 8(K)] Oly (gy | + UMAR UK) (4.10) 


§ _ (k) 


a 


The optimal conte lw.) eae 


u(k) = ~ [k, k, k, kK, |¢ = - K e(k) (4.11) 


where K is the steady-state gain matrix. The actual control 


u(k) is then obtained by integrating u(k) to obtain: 


E,, (m) 
k k pD(m) 
(k) = (m) = - (ke Key ee eee = 
u y u(m — ds 42 ead § , (m) (4.12) 
6, (m) 
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The LOR weighting matrices, Q and R, are chosen to satisfy 
the following design criteria: 


1. The overshoot to a step input should be less than five 
percent. 


2. The five percent settling time, t;,, is less than or 
equal to 1 second. 


Using an iterative process, the weighting matrices that meet 


the above design criteria were found to be: 


Osmo 20) 0 

0 000 

2 eae (4.13) 
Slo 0008 Re(1) 

0 a0 0 


The resulting steady-state feedback gain matrix, K, is: 
K=[k, k, k, k,} = [0.5347 -0.2385 0.1329 0.0017] . (4.14) 


Figure 4.1 shows the roll rate controller block diagram. 
d. Simulation Results 

Figure 4.2(a) shows the response of the closed 
loop system to a step input of six degrees/second. The 
design criteria of an overshoot less than five percent and 
the five percent settling time, t,,, less than one second 
are achieved. Figure 4.2(b) shows the aileron vane angle as 
a result of the step input. Weir (Ref. 16] demonstrates 
that the control vanes stall when displaced by an angle of 


plus or minus 30 degrees from the air flow zero reference. 


57 


Figure 4.1 Full Order Roll Rate 
Controller Block Diagram 
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Full Order Controller step Response 


Figure 4.2 


eTbuy eurwa (q) 
(soos) ouity 
sv Y st 3 cic z cl I 50 
2 Bae a 





a[suy ue, uolary 
esuodsey deqg (v) 
(soos) awit, 


Sb bv St t 4 c ST l $0 


b86'E= 100YSIMO % 





: Se ee ee ho —— 1S 


ey []Oy 





(Zap) spniqdury 


(d08s/8ap) apniydury 


59 


Therefore, it is important to ensure that the design not 
require vane angles greater than 30 degrees. 

Because the roll rate controller was designed 
using the linear model, it is necessary to measure the 
relative robustness of the design to ensure that it will be 
applicable to the nonlinear model. The phase and gain 
margins are such measures. From Figure 4.3, the gain margin 
1s equal to 32.63 dB and the phase margin is equal to 63.99 
degrees. A general rule of thumb is to require a gain 
margin greater than six dB and a phase margin greater than 
30 degrees to ensure satisfactory performance. Clearly, 
these benchmark values are exceeded and the above design 
should perform well when applied to the nonlinear model. 

The steady-state gains of Equation (4.14) are an 
optimal solution for the performance measure of Equation 
(4.10). However, this optimal solution requires that all of 
the states be measured or estimated. Although measurement 


of the vane angle, 54 is possible, a computational observer 


al 
would have to be included to provide an estimate of the vane 


velocity, 6,. Inclusion of an observer would add an 


additional task to the onboard digital computer, and 
increase the complexity of the controller design. An 
alternative to the observer would be simply to set the vane 
velocity gain to zero. The system would have to be 


Simulated with this gain equal to zero and the phase and 
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Full Order Controller 


Gain and Phase Margins 


Figure 4.3 
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gain margins computed to determine how the system response 
would be affected. 
e. Reduced Order Model 
Because the servo dynamics are significantly 
faster than all other system dynamics, a second alternative 
is to ignore the servo dynamics completely and form a 


reduced order state space representation of the system. 


Eee 
pt} |o 0 


The discrete-time state equation is given as 


ED 


BD 


Anal (4.15) 








0 
+ 
hell 











— / - -0.01 
poke + 1.) 0 mC, 





~~ eA: i). 
D(k) 0 .2izeiee- 





Now, the steady-state gain matrix, K, is determined from 


Equation (4.16) and the performance measure: 





z Erealphc 
wi = {tz ta p(k)]Q P| : +uT(kK)R auf _ (4-217) 
k=0 D(k) 
The optimal control wu(xkK) is: 
Pade 
Ui) "= [ oe = -Ke(k) . (4.18) 





and u(k) is 
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= =e ipim) | (4.19) 


The weighting matrices, Q and R, are defined as 


{0.3 0] v2 (4.20) 
of, peta. 


The resulting steady-state feedback gain matrix, K, is 
Feenent3 76 —02027) . (4.21) 


Figure 4.4(a) shows the response of the system to 
a step input of six degrees/second. The design criteria for 
overshoot and settling time are achieved. The step response 
for this reduced order model is almost identical to the 
response for the full order model given in Figure 4.2(a). 
From Figure 4.4(b), the maximum magnitude of the aileron 
vane angle displacement is 0.47 degrees. This vane 
displacement angle meets the requirement of less than or 
equal to 30 degrees. 

The necessity of determining the relative 
robustness of the system is now magnified due to the reduced 
order controller design. From Figure 4.5, the gain and 
phase margins are equal to 20.59 dB and 57.66 degrees 
respectively. Note that these margins were computed for the 
controller with the full order model. The reduction in gain 


and phase margins from the full order model is an expected 
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Reduced Order Controller 


Phase and Gain Margins 


Figure 4.5 
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qP urey 


consequence of the reduced order controller. However, the 
requirement for a minimum of six dB gain margin and 20 
degrees of ohase margin is satisfied. Therefore, it is 
expected that the reduced order controller design will be 
applicable to the nonlinear model. Figure 4.6 is a 
graphical realization of the reduced order =zoncrolles 
applied to the roll rate system. 

f. Summary 

An optimal tracker was designed for the roll rate 
control by applying the procedures of Chapter III. This 
control system yielded very satisfactory system performance 
and robustness. Recognizing the dynamics of the servos 
could possibly be ignored, a second reduced order design was 
developed. The resulting controller proved to be robust and 
almost identical in performance to the full order controller 
design. 

Because of the favorable results obtained with 
the reduced order roll rate controller design, the altitude 
rate and the pitch and yaw angle controllers will be 
designed using this technique. The resulting altitude rate 
and pitch and vaw controllers will be evaluated analogously 
to the roll rate controller to ensure they meet or exceed 


the design criteria. 
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Figure 4.6 Reduced Order Roll Rate 
Controller Block Diagram 
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C. ARCHYTAS ALTITUDE RATE CONTROLLER 
1. The Altitude System 

The primary function of the aititude rate Ccontromees 
for Archytas is to allow the operator to command a desired 
translational velocity along the vehicle’s longitudinal, oF 
x axis. This translational velocity wii Senile eae 
operator to position Archytas at a given altitude by 
terminating the velocity when a desired altitude is 
achieved. This capability is critical during the landing 
phase in order to control the altitude as Archytas lands in 
a vertical position. 

The simplified equation of motion which describes 


Archytas’ altitude rate subsystem is given as: 


i= h = 0. 10029507 (4.22) 


Note the change of notation from wu to 4, which is more 


appropriate. The throttle servo model is identical to the 


aileron servo model and is given as: 
6.=-74.16,.-2745.86,+2745 Sam (4.238) 

The Archytas engine is modeled as a first order ag: 

ent Olimen ope (4.24) 


The altitude rate tracking erromeis)det ined tas 
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2,2f =e) ; (4.25) 
where 4, is the input command and 4 is the measured 


altitude rate. From Equation (4.25), the differential 


equation for the tracking error is: 


B.=-A . (4.26) 


Sembining Equation (4.26) swith the time derivatives of 
Bquations (4.22), (4.23) and (4.24), a state equation that 


is appropriate for tracking system design is: 


“AY 19 =a 0 0 “1 0 
A 0 0 -0.10029 0 0 A 0 
eo 2) tees 65 OO + | Olt (4-27) 
5. 0 0 0 0 1 6. 0 
§ 10 0 0 -2745.8 -74.1)|/% | (2745.8) 
t | 


In a manner identical to that used for the roll rate 
controller design, the reduced order altitude rate model is 
formed. The controller is designed using the reduced order 
model and its performance and robustness are evaluated with 
respect to the full order model. Ignoring the servo 


dynamics, the reduced order state space model is: 


Ey Coa 0 Ey 0 
| = |\Ome0 0.10029) )| A | + 0 iio (4.28) 
ont 19 0 -2 Jfdzn) (2675.5 
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2. Altitude Rate Controller Design 
a. Discretizing the Altitude Rate System 
With the sampling period, T, equal to 0.01 


seconds, the discrete-time state equation is: 


Si el E(k) | = 
B(K+1) 0) = 10% gat) “Qqg0-0 tees ee (4.29) 
ba (Ktl) | 10 50 Oita ais emeeiela) Provan" 





b. Gain Determination 


altitude rate control loop are chosen to satisfy the 
following design criteria: 


1. The overshoot to a step input should be less than five 
percent. 


2. The five percent settling time, t;,, 1s less than or 
equal to 2 seconds. 


The resulting weighting matrices are: 


0) saa 0 0 
O=|} 0 0.02 0}; &={7] (4.30) 
0 0 0 
The steady-state feedback gain matrix, K, 1s: 
K=(-0.3060 0.2003 0.0038] - (4 ee) 


Figure 4.7 shows» the altitude rate cemtreller block diagram. 
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Figure 4.7 Reduced Order Altitude 
Rate Block Diagram 





Tak 


c. Simulation Results 

Figure 4.8(a) shows the response of the aititude 
rate closed loop system to a step input of ten feet/second. 
The design criteria of an overshoot less than five percent 
and the five percent settling time, t,., less than two 
seconds is achieved. Figure 4.8(b) shows the control 
applied to the system. 

The gain and phase margins were determined as a 
measure of the system robustness. From Figure 4.9, the gain 
margin is equal to 18.03 dB and the phase mardin is equal to 
58.78 degrees. These values of gain and phase margin 
indicate that the altitude rate control system design based 
on the reduced order model should provide good results when 


applied to the nonlinear system model. 


D. ARCHYTAS PITCH AND YAW ANGLE CONTROLLER 
1. The Pitch and Yaw Angle System 

The purpose of the pitch and yaw angle controller 
for Archytas is to maintain commanded orientation around the 
pitch and yaw axes. This requirement is necessary to allow 
the operator to position Archytas during landing pv pitching 
or yawing the vehicle slightly to induce 2 translatlenmam 
velocity along the body fixed y or z axis. This control 
problem is complicated by the gyroscopic coupling between 


the axes caused by the propeller. 
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Reduced Order Altitude Rate 


Controller Gain and Phase Margins 


Figure 4.9 
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The simplified equations of motion that describe the 


pitch and yaw angle subsystem are given 4s: 


fo ee oO Oar; (4.32) 


=-12.286,+1.78¢q . (4.33) 


Mahe elevator and rudder servos are modeled in a manner 
identical to the aileron servo and the corresponding 


differential equations are: 


Wee on 2745.00.+2745.8u, ; (4.34) 

and 
pea 7 47107-2745. 86.40945.8u, . (45:35) 
The pitch and yaw angle tracking errors are defined 

Ao. 
Ey=8.-8 ; (4.36) 

and 
Ts itr: ae (4.37) 


where 6. and Y. are the input commands and 9 and YP are the 


measured pitch and yaw angles, respectively. From Equations 
(4.36) and (4.37), the differential equations for 


the tracking errors are: 
Ey=-q | (4.38) 


and 


SD 


: -_ (4.39) 


Combining Equations (4.38) and (4.39) with the time 
derivatives of Equations (4.34) and (4.35), a state equation 
that 1s appropriate for tracking system design is X=A XtBY, 


where the state is: 


xT=([Ey Ey d Ey By t 6, 6, 6, 5, 17 (4.40) 


loam] 
= 


the control matrix is 1’=!u, u.])“, and the system matureae 


ace : 
(0) Chae ol 0 0 0 0 0 
0 <1 Gad 0 0 0 0 0 
00 0 0 0 Weaic2iem ae 0 0 0 
0 0) ee 0 0 0 0 0 
aero 0 Oa 0 0 0 0 (4.41) 
#10 01.78 0 0 0 0 ~12.28 0 0 
mo 0 “onc 0 0 0 1 0 
m0 Coe 0 0 0 0 S 
7 Oo OF 0G 0 -2745.8 0 -74,.18 0 
loo 0 00 0 0 -2745.8 0 “era. 
and 
oh 
pr (09000000 2745.8 0 | | (4.42) 


Ome 0.0 0 Oui 0 21a Se 


In a manner similar to that used in the roll sae 
and altitude rate controller designs, the reduced order 
model is formed by ignoring the servo dynamics. The 


resulting state is: 
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ees CO La Sa Fhe ; (4.43) 


mite, CONELOP Matrix is v-={6 





are: 
(eee. OO 
(me = 1 ~ Oman 0 
fGeeeO eie Oe =1.52 
FP = ; (4.44) 
— ro Oo ~ 1° 0 0 
. — 0 
e 6) = 0 aa) oO wl 
and 
fh 
Mcreo ins 00 860 | (4.45) 
Sen: 0 0 ee oe nen ia 


2. Pitch and Yaw Angle Controller Design 
a. Discretizing the Pitch and Yaw Angle System 
Using the sampling period, T, equal to 0.01 
seconds, the discrete-time state is: 
x(k) 7 = [Eg(k) £y(k) @(K) Ey(k) Ey(k) r(k)]"; (4-46) 
the control matrix is yv(k)7 = [6,(k) 4 


aided’ geaga the 


discrete-time system matrices are: 


TT 


Fe 0.01 9 0 ‘ 0 | 
10062) 0.62 0 0 ma O10) 
0 0 0.9999 0 0 -0.0162 
® (k) ' ie. | (4.47) 
ioumeuass ee 0 
0 0 0 0 0.99 0 | 
O» Oo CnONOZiGRaaO 0 0.9999 | 
and 
0 i | 
0.0006 0 | 
| 
-0.1119 0.001 
igen / (4.48) 
0 0 
S50 0 





| -G.001 -0.2228] 


b. Gain Determination 
The LQR weighting matrices, Q and R, for the 
pitch and yaw angle control loop are chosen to satisfy the 
following design criteria: 


1. The overshoot to a step response should be less than 50 
percent. 


2. The five percent settling time, t;,, is less than or 
equal to 2 secends. 


The allowable overshoot requirement may seem too liberal. 
However, due to the small pitch or vaw angles (one to five 
degrees) required to induce translational velocities along 
the body fixed y and z axes, the overshoot should not 
present any problems to the operator. Computer simulations 


have indicated that excessive control values, and vane angle 
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deflections would be required, if the controller was 
designed to limit the overshoot to values less than 15 
percent. The large overshoot values are due to the reduced 
Seeer COntroller. If a full order controller is 
implemented, the large overshoot values could be eliminated 
with small vane angles and a reasonable amount of control. 
However, that would require additional sensors and/or an 
Observer as discussed previously. 

The resulting state weighting matrices were found 


to be: 


6.0 0 ONCE 0 

= Gone G 
g=|9 9 9500 OF i ia (4.49) 
> io 0 6 om = Sona 

bao O80 1 0 

lo ip Ome C ost 


The steady-state feedback gain matrix, K, is: 


_| 2.2741 1.8933 -0.6059 0.6390 0.5320 we (4.50) 
Once Sn os20 -Oseg24 @.2741 1.8933 -0.6059) 


Figure 4.10 shows the MIMO pitch angle and yaw angle 
controller block diagram [Ref. 14: p. 179]. 
c. Simulation Results 
Figure 4.11(a) show the response of the system to 
a commanded pitch angle of ten degrees and a commanded yaw 
angle of zero. The large value of overshoot was expected 


due to the severe coupling of the axes. Figure 4.11(b) 
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Figure 4.10 Reduced Order Pitch Angie and 


Yaw Angle Controller Block Diagram 
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Figure 4.11 
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(s8ap) opriqdury 
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(b) Yaw Angle 


shows the induced yaw angle from the commanded pitch angle. 
Figure 4.12 gives a Similar representation. The yaw angle 
is commanded to five degrees while the pitch angle is held 
at zero. Finally, Figure 4.13 shows the response when the 
pitch and yaw angle are commanded simultaneously to five 
degrees. The maximum vane deflection applied in the 
Simulations documented by Figures 4.11, 4.12, and 4.13 was 
1.7 degrees. The design goals for overshoot and settling 
time are satisfied for the linear model with reasonable vane 
deflections. 

d. Singular Value Analysis 

Evaluating the robustness of the roll and 

altitude rate controllers was accomplished by computing 
their respective gain and phase margins. In the case of the 
MIMO pitch and yaw control system, the description of 
robustness in terms of the gain and phase margins becomes 
more complex. Therefore, a different approach is taken in 
order to determine the robustness of the pitch and yaw angle 
control system. The method of choice is to apply singular 
value analysis to the system with perturbations. 
Maciejowski (Ref. 19] presents a detailed development of the 
theory and procedures for computing the singular values of a 


MIMO system and using them for robustness analysis. 
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(b) Yaw Angle 
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Singular value analysis is used to examine the 
robustness of the MIMO pitch and yaw angle controller. Burl 
(Ref. 13] demonstrates that a system is stable for all input 


multiplicative perturbations, 4, such that 


(4.51) 


where o is the largest singular value of the matrix, and 


T(jw) is the transfer function from the output of the 
perturbation to the input to the perturbation as shown in 
Figure 4.14. The singular values for the MIMO pitch and yaw 
angle system are plotted in Figure 4.15. The largest 


Singular value, o, is equal to 1.5. Therefore, applying 


Equation (4.51), the MIMO system is expected to be robust 
for perturbations such that the infinity norm of 4 is less 
than or equal to 0.667. This indicates the controller 
design should provide good results when applied to the 


nonlinear model.[{Ref. 13] 


E. RESULTS WITH THE NONLINEAR SYSTEM 

The controllers designed in the previous sections were 
designed using optimal control theory. Because the servo 
states were ignored in determining the gains, the steady- 
state gains are sub-optimal in respect to the full order 
system models. Through simulations, and performance and 


robustness evaluations it was determined that ignoring the 
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servo dynamics would have no adverse affects. This section 
will apply the controllers designed with respect to the 
linear system models to the nonlinear model. The MATLAB 
programs used to simulate the nonlinear system are included 
as Appendix A. 
1. Simulation One - Figure 4.16 

The system was commanded to climb at a rate of ten 
feet/second for five seconds. The slope of the altitude 
plot indicates that the rate of ten feet/second is achieved. 
The roll rate, pitch angle and yaw angle were set equal to 
zero. The control system is tracking the commanded altitude 
rate input while regulatingsthe other inputs to 7eree 
Figure 4.16 shows that the vehicle climbs to an altitude of 
50 feet. Note the coupling between the roll axis (both roll 
rate and roll angle) as the engine accelerates. The 
displacement of the aileron control vanes do not exceed plus 
or minus one degree in deflection. 

2. Simulation Two - Figure 4.17 

Figure 4.17 shows the response to a commanded 
altitude rate of 25 feet/second for four seconds. A 
commanded roll rate of ten degrees/second for two seconds is 
applied starting at time equal to six seconds. The slope of 
the altitude plot indicates that the altitude rate of 25 
feet/second is achieved. The vehicle overshoots the desired 


100 feet altitude and begins to approach a steady-state 
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value of approximately 80 feet. The large overshoot for the 
altitude is due to the momentum of the vehicle. Climbing at 
a high rate of climb, the vehicle is expected to exhibit 
this type of behavior. Additionally, the control system is 
controlling the altitude rate, not the altitude. The 
operator would have to account for this in his control 
adjustments during flight. Again, the coupling between the 
roll axis and the engine acceleration is observed. The 
commanded roll rate drives the vehicle to an roll angle of 
18°. The two degree error is due to the minus two degrees 
of roll angle caused by the reactive torque to the altitude 
rate. The increased aileron vane displacement in order to 
perform the desired maneuver is shown. | 
3. Simulation Three - Figure 4.18 

The vehicle was commanded to climb at a rate of 50 
feet/second for four seconds. A pitch angle of minus five 
degrees was commanded at time equal to six seconds. Figure 
4.18 illustrates that the increased rate of climb increases 
the altitude overshoot. The coupling between the engine and 
roll axis is present. Additionally, the coupling between 
the pitch and yaw axis is shown. As the vehicle is pitched, 
the yaw angle is perturbed. Note, the controller is 
limiting the eee between the pitch and yaw axis. The 
yaw angle is perturbed less than one degree for a commanded 


pitch angle of minus five degrees. 
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Figure 4.18 Nonlinear Simulation Three 
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4. Simulation Four - Figure 4.19 
Figure 4.19 shows the system response for an 

altitude rate of 50 feet/second for four seconds. The 
commanded roll rate is ten degrees/second for one second at 
time equal to three seconds. The commanded pitch angle is 
minus three degrees applied at time equal to six seconds. 
The control vanes displacement angles all remain less than 
four degrees for this series of maneuvers. The controller 
is able to provide satisfactory results for the nonlinear 


model. 


F. CONCLUSION 

The controller has demonstrated the ability to control 
the nonlinear system with varying inputs applied. It can be 
concluded that the robustness of each control system is 
sufficient in order to account for the difference in system 
parameters between the linear models and the nonlinear model 
of Archytas. The fact that the control systems maintain 
stability while directing commanded inputs to steer and 
maneuver the vehicle supports the linear modeling approach 
taken in this thesis. The next step in the design process 
is to test the control systems on a prototype vehicle in a 
controlled experiment. Such a test was performed with the 
roll rate control system and an Archytas prototype mounted 
on a test stand. The results of this test will be discussed 


in Chapter V. 
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V. CONCLUSIONS 


A. ROLL RATE CONTROL SYSTEM FIELD TEST 

Figure 5.1 shows the Archytas prototype (minus the wings 
and canard) mounted on a test stand that allows the vehicle 
to spin about the longitudinal, or x axis. This 
configuration was used during the testing of the roll rate 
controller. The testing was accomplished in a qualitative 
manner; no empirical data was collected, the object of the 
test being to validate the roll rate controller design. The 
test consisted of three parts: 

1. The ability of the roll rate controller to eliminate 
the rotational velocity imparted to Archytas from the 
reactive torques applied to the roll axis as the engine 
speed is varied (decoupling). 

2. The ability of the roll rate controller to eliminate 
the rotational velocity imparted to Archytas from the 
effect of cross-winds (disturbance rejection). 

3. The ability of the roll rate controller to allow the 


operator to impart or terminate a rotational velocity 
about the x axis (tracking commanded inputs). 


Part one was accomplished by varying the engine speed 
from idle to maximum rpm with the commanded roll rate equal 
to zero. Each time the engine speed was varied, the 
rotational velocity of the vehicle was eliminated by the 
roll rate controller. It was observed that the rotational 


velocity was eliminated within approximately two seconds. 
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Archytas Prototype Mounted 
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on the Test Stand 


Figure 5.1 


This response is consistent with both the linear and 
nonlinear Simulations. 

Part two was accomplished by setting the engine speed at 
constant values ranging from idle to maximum rpm with the 
commanded roll rate equal to zero. A disturbance about the 
X axis was imparted to Archytas by pushing the test stand 
mounting bracket connected to the vehicle. The effect of 
pushing the connecting bracket is analogous to a rotational 
velocity imparted by cross-winds. It was observed that the 
rotational velocity was eliminated within approximately two 
seconds. This response is consistent with both the linear 
and nonlinear simulations. 

Part three was accomplished by inputting a commanded 
roll rate while the engine speed was held constant and also 
varied. It was observed that the operator could position 
Archytas at a desired angle under either test condition. 
Additionally, disturbances were imparted to Archytas during 
these tests. The roll rate controller proved robust enough 
to eliminate the disturbances and allow the operator to 
position the vehicle. This response is consistent with both 
the linear and nonlinear simulations. 

Based on the computer simulations, both linear and 
nonlinear, and the results of the above tests, it is 


concluded that the roll rate controller design is valid. 
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B. FUTURE RESEARCH 

The development of the Archytas control systems was 
accomplished using the data from the AROD as first 
approximations. Future research will include the 
computation of Archytas specific data through empirical 
methods or obtaining the data by wind tunnel testing. This 
would allow for a better model of the Archytas to be 
developed. 

The second order dynamical model of the servos was 
obtained with the vanes mounted on the servos with zero 
downwash from the propeller. Future research should include 
modeling of the servos with the engine at hover rpm and the 
vanes positioned within the downwash. 

The effects of inclusion of the vane position angles in 
the control systems could be evaluated further. [In 
particular, the effect the vane angle position feedback 
would have in the reduction of the overshoot in the pitch 


and yaw angle controller should be investigated. 
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APPENDIX A 


MATLAB SIMULATION PROGRAMS 


Ooeooa$acwaoecegcocncoeoocoeooeoeooeoeoeoo oo ooeaeaetoooeceoaoceoeceeoooeooooooO oOo oOo ood 8 
BCCCCCCCCRCCCEREDEREEEECERECREERERERRERREEEEECCEEECEEEEEEEESEEEEES 
PROGRAM NAME nonlin.m 
e 
DESCRIPTION: 


This program will simulate the Archytas nonlinear model 
with the linear controllers. The Control Laws are the 
Steady-state Linear Quadratic Regulator solutions with 
the necessary modifications to the performance measure. 
CECEEECCEECCCCCEEECCECEEECEECEEEESEEEEEEEEEEEESEEECECEEESS 
ar 
Call the program noninp.m to enter the desired roll rate, 
alttitude rate, pitch and yaw angles. 
oninp 
Call the program in _cond.m to load the initial conditions 


oe J © oP oP oO oP oP oN? 


° 
© 


oi? 


2 
6 
e 


Of TF oP OP TQ oP oP oP oP DO oP oP oO oO 


and * constants for the simulation. 

in_cond 

* Call the program n_1]_inp.m to load the wind tunnel data. 
nel inp 


% Loop for simulation 
for k=l:num_of stepsl 


* Compute the speed of the propeller 
Speed(k) = speed hover + del rpm(k); 
%* Compute the thrust supplied by the engine to the vehicle 
thrust(k) = thrust_hover + Xrpm * del rpm(k) ; 
6 Check thrust limits 
if thrust(k) < thrust min 
thrust(k) = thrust min; 
elseif thrust(k) > thrust_max 
Euruse(k) = thrust.max; 
end 


velocity tip(k) = speed hover + del rpm(k); 
del tip(k) = velocity tip(k) - speed hover; 
tip _del(k) = velocity tip(k) - speed(k); 

stad (k) = tip del(k); 


% Compute vehicle velocity and angle of attack total 
aoatot(k) = alpha max; 
velocity tot(k)=sqrt(u(k) *2+v(k) *2+w(k)%*2); 
if velocity tot(k) > vaero 
vwterm = sqrt(v(k)*2 + w(K)%*2); 
aoatot(k) = asin(vwterm/velocity tot(k)); 
if aoatot(k) < alpha min 
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aoatot(k) = alpha min; 
elseif aoatot(k) > alpha max 
aoatot(k) = alpha max; 
end — 
end 


veq(k) = tablei(veq_table,aoatot(k))j; 

newveq(kK) = sqrt(weight ratio) *veq(k) ; 

rhova(k) = rho * newveq(k) * pi; 
velocity tip 2(k) = 2.0/velocity tip(k); 
velocity delta(k) = velocity tot(k) - newveq(k); 


% This section assumes that no aerodynamic forces and moments 
$ exist 
if velocity tot(k) <= vaero 


fax 

fay 

faz 

mrx 

mpy 

my Z 

patch jeeam 
yaw_ trim 
prech tacter 
yaw _factr 


hrust(k) * gravity/weight; 
.0; 


eo ™e@ 


t 
0 
0 
0 
0. 
0 
0 
0 
1 
Ji 


Ooo oo0o00 Oo 


elseif velocity tot(k) > vaero 
if aoatot(k) <= alpha _min 


thrust(k) * gravity/weight; 
0.0; 


fax 

fay 

faz 

mrx 

Mpy 

my Z 
PLtchetEan 
yaw _ trim 
Peechweracer 
yaw factr 

end 


—- © © oe Oro 
Oooo oo00 0 
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else 


% This section computes the aerodynamic forces and moments 


lleq = weight ratio * leq; 
ddeq = weight ratio * deq; 
sseq = weight ratio * seqq; 
req = tablel(req_table,aoatot(k)); 
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peq ~suamlol(peq table,acatot(k) ) ; 
yeq = Gablhel(yeq table,aoatot (k).) 
rreq = weight ratio * req; 

ppeq = weigth ratio * peq; 

yyeq = weight _ ratio * yeq; 

eel stan le l(eldelwcable,acatot(k),) ; 
cddel = tablel(cddel table,aoatot(k)); 
csdel = tablel(csdel table,aoatot(k) ) ; 


ldel = velocity tip 2(k)*lleq-rhova*weight ratio*cldel; 
ddel=velocity tip 2(k) *ddeq-rhova*weight ratio*cddel; 
sdel =velocity tip 2(k)*sseq-rhova*weight ratio*csdel; 


crdel = tablel(crdel table,aoatot(k) ); 
cpdel = tablel(cpdel table,aoatot(k) ); 
cydel = tablel(cydel table,aoatot(k)); 


rdel=velocity tip 2(k)*rreq-rhova*weight ratio*crdel; 
pdel =velocity tip 2(k)*ppeq-rhova*weight ratio*cpdel; 
ydel=velocity tip 2(k) *yyeq-rhova*weight_ratio*cydel; 


lslope = tablel(lslope table,aoatot(k) ) ; 

dslope = tablel(dslope table,aoatot(k)); 

sslope = tablel(sslope table,aoatot(k) ); 

SE ] = lleqtidel*del tip(k)+lslope*weight_ratio* 
velocity delta(k); 

fia = ddeq+ddel*del tip(k)+dslope*weight_ratio* 
velocity delta(k) ; 

ffs = sseqtsdel*del tip(k)+sslope*weight_ratio* 


velocity delta(k); 


= Jeaplel(Estopertapre,qagaror (kK) ); 
pslope = tablel(pslope table,aoatot(k)); 
= yi 


yslope Eapleinvetopegtable ,aoatot (k 

mr = rreqtrdel*del tip(k)+rslope*weight_ratio* 
velocity delta(k); 

mp = ppeqt+tpdel*del tip(k)+pslope*weight_ratio* 
velocity delta(k) ; 

my = yyeqtydel*del tip(k)+yslope*weight_ratio* 
Velocity deliask) ; 

fs = ffs * gravity/weight; 

fia: = ffl * gravity/weight; 

fd = ffd * gravity/weight; 

delta = atan(v(k)/w(k)); 

fax = fl * sin(aoatot(k)) - fd * cos(aoatot(k) ); 
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fayz =-fl * cos(aoatot(k)) - fd * sin(aoatot(k)); 
fay = fayz * sin(delta); 
faz = fayz * cos(delta); 
mrx = my * sSin(aoatot(k)) - mr * cos(aoatot(k)); 
mpyyz =-my * cos(aoatot(k)) - mr * sin(aoatot(k) ); 
mpy = mpyyz * sin(delta); 
myzZ = mpyyz * cos(delta) ; 


velocity trim =tablel(velocity trim _table,aoatot(k1)); 
vaneff = tablel(vaneff table,aoatot(k) ); 


pitch trim = velocity trim * cos(delta); 
yaw trim =-velocity trim * sin(delta) ; 
pitch factr = vaneff; 
yaw factr = vaneff; 
end 
lat(k) = ((rscale*mrx*gravity*144.0) /Ixx)+(La*del_a(k))- 


(velocity delta(k)*2 /Ixx) * prop _torqd; 


mat(k) = ((mpy*gravity*144.0) /Iyy)+(Me*pitch factr* 
(del _e(k)-pitch_trim)); 
((myz*gravity*144.0)/1I1zz)+(Nr*yaw factr * 
(del _r(k)-yaw_trim)); 


nat (k) 


° 


% Begin Differential Equations of Motion 


% Altitude Differential Equations 

temp1 (k) =u(k) *cos(theta(k) ) *cos(psi(k) ); 

temp2 (k)=v(k) *(sin(phi(k) ) *sin(theta(k) ) *cos(psi(k))- 
cos (phi(k))*sin(psi(k) )) ; 

temp3 (k) =w(k) * (cos (phi(k) ) *sin(theta(k) ) *cos(psi(k))+ 
Sin(phi(k))*sin(psi(k))); 

u_earth(k)= temp1(k) +ttemp2 (k) +temp3 (k) ; 


% Heading Angle Differential Equations 

temp5 (k)=p(k) *cos(theta(k) ) *cos(psi(k)); 

temp6 (kK) =q(k) *(sin(phi(k) ) *sin(theta(k) ) *cos(psi(k))- 
cos (phi(k) )*sin(psi(k))); 

temp7 (k) =r(k) * (cos (phi(k) ) *sin(theta(k) ) *cos (psi (k) ) 
+Ssin(phi(k) )*sin(psi(k))); 

hdg_ dot (k)=temp5 (k) +temp6 (k) +temp7 (k) ; 


% Euler Angle Differential Equations: Pitch, Yaw, and Roll 

theta dot(k) =q(k) *cos(phi(k) )-r(k) *sin(phi(k) ); 

phi _dot=p(k)+(q(k) *sin(phi(k) )+r(k) *cos(psi(k)))* 
tan(theta(k)); 
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psi_dot(k)=(q(k) *sin(phi(k) )+r(k) *cos(phi(k)))/ 
cos(theta(k)); 


% Vehicle Pitch, Yaw and Roll Rate Differential Equations WRT 
% Body Axes 

p dot(k) = (bca*q(k) *r(k))-(Iralph*stad(k))+lat(k) ; 

(cab*p(k) *r(k) )-(Irx*r(k) *speed(k) /Iyy) +mat (k) ; 
(abc*p(k) *q(k) )+(Irx*q(k) *speed(k) /1zz) +nat(k) ; 
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% Vehicle Velocity Differential Equations WRT Body Axes 

u_dot (k) =( (v(k) *r(k) )-(w(k) *q(k) )) -gravity*cos (theta(k) )* 
cos (psi(k) )+fax; 

v_dot (k) =( (w(k) *p(k) )-(u(k) *r(k) ))+gravity*cos (theta (k) ) * 
sSin(psi(k))+fay; 

pmoee( kK) — (i(k) eek) ) -(v(k)*p(k)))-gravity* 
Sin(theta(k))+faz; 


Define Error Variable in Pitch and Yaw 
_theta(k)=theta_com(k) -(theta(k)+noise(k) ) ; 
_psi(k)=psi_com(k)-(psi(k) +noise(k) ); 


tH by ol 


Define Error Variable in Roll 
p(k) =p_com(k) -(p(k)+noise(k) ) ; 


fr oe 


Define Error Variable in Altitude 
_h_dot(k)=h_dot_com(k) -(u(k)+noise(k) ) ; 


ty of 


Integrate Errors 
@eketain(k+1l)—E wehetain(k)+Ts*E theta(k) ; 
_psiin(k+1)=E_ psiin(k)+Ts*E_psi(k) ; 
_pin(k+1)=E_pin(k)+Ts*E_ p(k); 
Pamdotingk+1)=E h dotin(k)+Ts*E h dot(k); 


Mmm wy oe 


Define Control System Gains 
Roll Rate Control Gains 
(eee 6 -O.2027/7) > 


% Pitch and Yaw Angle Control Gains 

fh —|| 222741 158933 —-0.6059 0.6390 0.5820 90.0024, 
meebo 90 —O85320 -0.0024 2.2741 1.8933 -0.6059,),; 

io—(-vesvo0 Of2003 0.0038]; 


% Calculate the Aileron (roll) Servo Control Input 
Ua (kK) =-K1*(E_pin(k+1) ; (p(k) +noise(k))]; 


Calculate the Elevator (pitch) Servo Control Input 


Ue (k)=-K2(1,:)*(E_thetain(k+1) ;E theta(k) ;q(k) ;E_psiin(k+1); 
E_psi(k) ;r(k)]; 
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¢ Calculate the Rudder (yaw) Servo Control Input 
Ur (kK) =-K2(2,:)*(E_thetain(k+1) ;E theta (k) ;q(k) ;E_ps1in(k+t1); 
E_ps1(k);r(k)]; 


6 Calculate the Throttle (altitude) Servo Control Input 
Ut (kK) =-K3*(E_h_ dotin(k+1) ;u(k) ;del_rpm(k) ); 


* Begin Integration of Equations of Motion 

del _a dot dot(k)=-H1*del a dot(k)-H2*del_a(k) +H2*Ua(k) ; 
del_e dot dot(k)=-H1*del_e dot(k)-H2*del_ e(k)+H2*Ue(k) ; 
del r dot _dot(k)=-Hl*del r dot(k)-H2*del r(k)+H2*Ur(k) ; 
del t dot _dot(k)=-H1l*del t _dot (k) -H2*del t(k)+H2*Ut(k) ; 


Roll, pitch, yaw 

p(k+1)=p(k)+Ts*p_ dot(k); 
qd (kil )=qik) +Tfs*qideotgay 
r(k+1)=r(k)+Ts*r_ dod; 


Velocities 

u(kK+1)=u(k)+Ts*u_dot(k) ; 
vt kt 1) =v (K) Esty dort) ; 
Ww (k+1) =w(k)+Ts*w_dot(k) ; 


éEuler angles 

theta (k+1)=theta (k) +Ts*theta” dot (k) ; 
phi(k+1)=phi(k)+Ts*phi _dot(k) ; 
psi(k+1)=psi(k)+Ts*psi_ dot(k) ; 


Servos 
del_a_dot(k+1)=del a dot(k)+Ts*del a dot dot(k); 
ipwaet a SaGiE (rs) > rmax 
dela _dot(kK+1) = rmax; 
elseif del a ~A6C(Kk+1) <=rirax 
del_a dot(k+1) = -rmax; 
end 


del _e dot(k+1)=del e dot(k)+Ts*del e dot dot(k); 
if del _@ dot(k+1) > rmax 


del e dot(kK+1) = rmax; 
elseif dele _Got(k+1) < -rmax 

Ged: e _dot (K+1) = -rmax; 
end 


del _r dot(k+1)=del r_ dot(k)+Ts*del r_ dot dot(k); 
if adem _Y doti@k+1) > Emax 


del r “GOt (kPl)"= max 
elseif del r _dot(k+1) < -rmax 

del r _dot (k+1) = —-Ymax; 
end 
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del t_dot(k+1)=del t_dot(k)+Ts*del t dot _dot(k); 
rf @e® t Beottukd 1) ) > =rmaxt 
del ie Beot( k+l) = rmaxt; 
elseif del t dot(k+1) < -rmaxt 
fel etesaot(k+t1) = -rmaxt; 
end 


% Engine 
del rpm dot(k+1)=-w_e*del rpm(k)+Ke*w_e*del t(k); 


Seed (k+tlj=del a(k)+Ts*del_ a dot(k) ; 


del e(kt+1)=del e(k)+Ts*del e dot(k); 
me dew Reger) > waxeftl 
del e(k+1) = maxdfl; 
elseif del _e(k+t1l) < -maxdf 1 
a@el e(k+tl) = -maxdfl; 
end 


del r(k+1)=del r(k)+Ts*del r dot(k); 
if del r(k+1) > maxdfl 
del r(k+1) = maxdfl; 
elseif del _r @@r1) < -maxdafl 
del _x(K+t1) = -maxdf1; 
end 


del t(kt+1)=del t(k)+Ts*del_t_ dot(k); 
if del me (k+l) >*tameax 
del t(k+1) = tamax; 
elseif del t(kt+1) < -tamax 
gel) t(kri) = -tamax; 
end 


% Engine 
del rpm(k+1)=del rpm(k)+Ts*del_ rpm_dot(k); 


%* Altitude 
gee (h+1l)—alct(k)+Ts*u_ earth (k) ; 


% Heading Angle 
heading (k+1) =heading(k)+Ts*hdg_dot(k) ; 


€(k)=Ts*(k-1) ; 
* End of simulation loop 
end 


€(kK+1)=Ts*k; 
u_earth(k+1)=u_earth(k) ; 
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% Plot the Results of the Simulation 

subplot (22 l=, paeot @teralt) ,@rid;xlabel( ta 
(sec)’),ylabel(’Altitude (ft) ’); 

title(’Altitude): 

subplot(222) ,plot(t,p*57.29958) ,grid,xlabel(" fa 
(sec)’),ylabel(’Roll Rate (deg/sec)’); 

title(’Roll Rate - p’); 


subplot(223) ,plot(t,phi*57.2958) ,grid, xlabei( 2a 
(sec)’),ylabel(’Angle (deg) ’); 

title(’Roll Angle Phi’); 
subplot(224),plot(t,theta*57.2958),grid,xlabel(’Time 
(sec)’),ylabel(’Angle (deg) ’); 

title(’Pitch Angle Theta’); 

meta nonl 

subplot(111); 


subplot(221),plot(t,psi*57.2958),grid,xlabel(’Time 
(sec)’),ylabel(’Angle (deg) ’); 

title(’Yaw Angle Psi’); 
subplot(222),plot(t,del a*57.2958) ,grid ,xYabe( Siam 
(sec)’),ylabel(’Angle (deg) ’); 

title(’Aileron Vane Angle’); 


subplot(223),plot(t,del r*57.2958)%qridy, xlabel @@miae 
(sec)’),ylabel(’Angle (deg) ’); 

title(’Rudder Vane Angle’); 

subplot(224) ,plot(t,del e*57.2958) ,grid, xlabel (sea 
(sec)’),ylabel(’Angle (deg) ’); 

title(’Elevator Vane Angle’); 

meta non2 

subplot ( 
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This program prompts the user for the desired 
Simulation inputs. 
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DESCRIPTION: % 
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% Prompt the user for the desired simulation time. 
disp(’ Enter the simulation time in seconds.’); 
Simulation time = input(’>>>> ’); 


% Prompt the user to enter to desired sampling time - Ts. 
disp(’ Enter the desired sampling time -~- Ts.’); 
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Ts.= input’ >>>> ’): 


% Compute the number of simulation steps based on the 
desired simulation time. 
num of stepsl = round(simulation time/Ts) ; 


% Prompt the user for the length of the altitude rate step 
input time. 

disp(’ Enter the length of the altitude rate input in 
seconds.’); 

step input timel = input(’>>>> ’); 


*¢ If step input entered is greater than simulation time 
entered, set them equal. 
if step input _timel > simulation time 
step input_timel = simulation time; 
end 


if step _input_timel ~= 0 
% Prompt the user for the desired magnitude of the step 
6 input in # of deg/sec. 
disp(’ Enter the magnitude of the altitude rate input in 
feet/sec.’) 
h_ dot = input(’>>>> ’); 
elseif step input timel == 0 
h dot = 0.0; 
end 


* Compute the number of step input steps based on the 
% desired step input time. 
num of steps2 = step input _timel/Ts; 

= round(num_of_ steps2); 


* Compute the number of zeros to be padded to the input 
* based on the length simulation time and length of the 
% desired step input time. 


num_of zerosl = num_of_stepsl - num _ of_steps2; 


* Define the commanded input vector 
h_dot_com = (h dot * ones(1,num_of steps2) 
zeros(1,num_of zerosl1)]; 


% Prompt the user for whether he wants to include a roll 
rate in the simulation. 
ans = input(’Do you want to include a roll rate in the 
Simulation? y/n [{y]: ’,’s/’); 
if isempty (ans) 
ans =/y" > 
end 
if ans ==’y’ 
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% Prompt the user for the start time for the roll rate 
step input. 

disp(’ Enter the start time for the roll rate input.’); 

start timel = input(’>>>> /); 


% Prompt the user for the length of the roll rate input. 

disp(’ Enter the length of the roll rate imput in 
seconds.’); 

length_roll rate = input(’>>>> /’); 


Compute the number of zeros to pad front of the 
command vector. 

num_of zeros2 = start _timel/Ts; 

num_of zeros2 = round(num of zeros2); 


oP oe 


% Computer the number of step input steps and number of 
% zeros to pad back of the command vector 


num_of steps3 
num of steps3 
num of zeros3 


lengemaroll rare ae, 
round(num of steps3),; 

num of stepsl- 

round((start timel+length roll rate) /Ts) ; 


* Prompt the user for the desired magnitude of the roll 

%* rate input in # deg/se. 

disp(’Enter the magnitude of the roll rate in 
degrees/sec.’); 

pl = input(’>>>> ’); 

% Define the commanded input vector 

p_com = [(zeros(1,num_of zeros2) pl * 0.017453 * 

ones(1,num_of steps3) zeros(1,num_of_zeros3)]; 


p_com = [(zeros(1,num_of stepsl)]; 
end 


% Prompt the user for whether he wants to include a pitch 
* angle in the simulation. 
ans = input(’Do you want to include a pitch angle in the 
Simtilativen?: y/i? jwise’ , so 
if isempty (ans) 
ans = ‘y’; 
end 
if ans ==’y’ 
% Prompt the user for the start time for the pitch angle 
% command. 
disp(’ Enter the time at which you desire to input the 
pitch angle.’) 
start time2 = input(’>>>> ’); 
* Compute the number of zeros to pad front of the 
% command vector 
num of zeros4 = round(start_time2/Ts) ; 
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Compute the number of ones to multiply the desired 
angle by to create the command vector. 

umber of onesl = num of stepsl - round(start_time2/Ts) ; 
Prompt the user for the desired pitch angle in 
degrees. 

disp(’ Enter the deSired pitch angle in degrees.’); 

theta com = input(’>>>> ’); 

theta _ com = ({zeros(1,num_of_ zeros4) theta_com * 0.017453 

ones(i,;number of onesl) ); 


of of FX ae ol? 


* 


else 
eheta com 
end 


([zeros(l1,num_of.stepsl) J; 


%* Prompt the user for whether he wants to include a yaw 
% angle in the simulation. 
ans = input(’Do you want to include a yaw angle in the 
Simulation? y/n [yj]: ’,’s’); 
if isempty (ans) 
ans = ‘vy’; 


if ans ==’y’ 
% Prompt the user for the start time for the yaw angle 
% command. 
disp(’ Enter the time at which you desire to input the 
yaw angle.’) 
start time2 = input(’>>>> ’); 
Compute the number of zeros to pad front of the 
command vector 
um of zeros4 = round(start_time2/Ts) ; 
Compute the number of ones to multiply the desired 
angle by to create the command vector. 
umber of onesl = num_of stepsl - round(start_time2/Ts) ; 
Prompt the user for the desired pitch angle in 
degrees. 
disp(’ Enter the desired yaw angle in degrees.’)j; 
DStgeeme= slinputd(/ >>ze—/.); 
psi_com = ({zeros(1,num_ of zeros4) psi_com * 0.017453 * 
ones(1,number of onesl) ]; 
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else 
psi _com 
end 


{zeros(1,num_ of stepsl) ]; 


Prompt the user for whether he wants to add sensor noise 
to the simulation. 


oP? oP 


ans = input(’Do you want to include sensor noise in the 
Similagion: y/n {y): ’,’s’); 
if isempty (ans) 
ans =/y’; 
end 
if ans ==/y’ 
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Define the random vector to represent the measurement 
noise. 

Assume the accuracy is +/- 1.745e-04. 
rand(’uniform’ ) ; 

noise = rand(1,num_of steps1l); 
vard = (1/12) * (1.745e-04)%*2; 
A sqrt (vard/cov(noise) ); 


of oP oP 


Desired variance 
Scalar multiplier to 
change the variance 


of oP oP 


noise = A .* noise; 
° 


% Adjust the mean to zero 


noise = noise - mean(noise) ; 
else 

noise = [zeros(1,num_of_ steps1)]; 
end 


% DESCRIPTION: 
% constant values used in no 


° 
% 
% 
This program inputs the initial conditions and % 
% 
% 
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* CONSTANTS 

vaero=0.5; speed max=859.0; speed min=649.0; 

prop _torq=0.0729; thrust_min=35.0; thrust_max=105.0; 
tamax=100.0; rmaxt=100.0; 

maxdf1=0.5236; rmax=0.87266; 

La=-21.04; Me=-9.01; Nr=-11.40; 

H1=71.1; H2=2745.8; 

w_e=2.0; Ke=837.758; Xrpm=0.2122; speed hover=712.0943; 
alpha_min=0.174533; alpha_max=1.570796; 

gravity=32.174; pi=3.1415962; rad_to_deg=180.0/p1i; 
Ixx=6908.4; Iyy=22944.64; I2z2=21515.64; Irx=41.62; 
bca=0.20685; cab=0.63663; abc=-0.74533; Iralph=0.00916; 
weight=85.0; thrust hover=85.0; velocity tip=722.5663; 
rho=0.00192; weight ratio=1.0; 

leq=85.0; deq=0.0; seqq=0.0; sslope=0.0801; rscale=0.0; 


INITIAL CONDITIONS 
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APPENDIX B CONTROL SERVOS 
The control servos incorporated into Archytas are 
position-commanded, quarter-scale model airplane servos. 
Five of these servos are used on Archytas: four to drive the 
aerodynamic control surfaces and one to drive the engine 
throttle. These servos can be modeled as second-order 


dynamical systems: 


ae 2 (B.1) 


R(s) s*+20w,+w% 





where C(s) is the output and R(s) is the input. ¢ and a, 
are referred to as the damping ratio and the natural 
frequency, respectively. 

The values for ¢ and w, were determined by experiment. 
The control vanes were mounted on the servos and position 
data was obtained for a range of step input commands. The 
input commands ranged from three to thirty degrees of vane 
deflection. By plotting the servo angle response data 
versus time, a series of step input response curves was 
developed. Using MATLAB, trial and error was used to fit a 
second order prototype system to the step response curves of 
the servos. Figure B.1 shows the "best fit" model. This 


model is defined by [=0.707 and w,=52.4 radians/second. 
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Thus, the second order servo model used in the design of the 


Archytas control systems is: 


Gis) 2745.8 (B.2) 


R(s) s°+74.1sS+2745.8 
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